      SUBROUTINE M28Y12(N, NZ, A, LICN, IRN, LIRN, ICN, U, IKEEP,IW, W,
     * LBLOCK, GROW, ABORT, IDISP, IFAIL, THR)
C     DERIVED FROM HARWELL LIBRARY ROUTINE MA28A
C     THROUGH F01BRF MARK 7 RELEASE. NAG 1978
C     AND Y12M -- ZLATEV & 0STERBY
C
C     THE PARAMETERS ARE AS FOLLOWS.....
C     N     INTEGER  ORDER OF MATRIX  NOT ALTERED BY SUBROUTINE.
C     NZ    INTEGER  NUMBER OF NON-ZEROS IN INPUT MATRIX  NOT
C     .     ALTERED BY SUBROUTINE.
C     A     REAL ARRAY  LENGTH LICN.  HOLDS NON-ZEROS OF
C     .     MATRIX ON ENTRY AND NON-ZEROS OF FACTORS ON EXIT.
C     .     REORDERED BY F01BRZ AND F01BRY AND ALTERED BY M30Y12.
C     LICN  INTEGER  LENGTH OF ARRAYS A AND ICN.  NOT ALTERED BY
C     .     SUBROUTINE.
C     IRN   INTEGER ARRAY  LENGTH LIRN.  HOLDS ROW INDICES ON INPUT
C     .     AND USED AS WORKSPACE BY M30Y12 TO HOLD COLUMN
C     .     ORIENTATION OF MATRIX.
C     LIRN  INTEGER  LENGTH OF ARRAY IRN.
C     ICN   INTEGER ARRAY  LENGTH LICN.  HOLDS COLUMN INDICES ON
C     .     ENTRY AND COLUMN INDICES OF DECOMPOSED MATRIX ON EXIT.
C     .     REORDERED BY F01BRZ AND F01BRY AND ALTERED BY M30Y12.
C     U     REAL VARIABLE  SET BY USER TO CONTROL BIAS TOWARDS
C     .     NUMERIC OR SPARSITY PIVOTING.  U=1.0 GIVES PARTIAL
C     .     PIVOTING WHILE U=0. DOES NOT CHECK MULTIPLIERS AT ALL.
C     .     VALUES OF U GREATER THAN ONE ARE TREATED AS ONE WHILE
C     .     NEGATIVE VALUES ARE TREATED AS ZERO.  NOT ALTERED BY
C     .     SUBROUTINE.
C     IKEEP  INTEGER ARRAY  LENGTH 5*N  USED AS WORKSPACE BY M28Y12
C     .     (SEE LATER COMMENTS).  IT IS NOT REQUIRED TO BE SET ON
C     .     ENTRY AND, ON EXIT, IT CONTAINS INFORMATION ABOUT THE
C     .     DECOMPOSITION.  IT SHOULD BE PRESERVED BETWEEN THIS CALL
C     .     AND SUBSEQUENT CALLS TO M28BSF OR F04AXF.
C     IKEEP(I,1),I=1,N  HOLDS THE TOTAL LENGTH OF THE PART OF ROW I
C     .     IN THE DIAGONAL BLOCK.
C     ROW IKEEP(I,2),I=1,N  OF THE INPUT MATRIX IS THE ITH ROW IN
C     .     PIVOT ORDER.
C     COLUMN IKEEP(I,3),I=1,N  OF THE INPUT MATRIX IS THE ITH
C     .     COLUMN IN PIVOT ORDER.
C     IKEEP(I,4),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN
C     .     THE L PART OF THE L/U DECOMPOSITION.
C     IKEEP(I,5),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN
C     .     THE OFF-DIAGONAL BLOCKS.  IF THERE IS ONLY ONE DIAGONAL
C     .     BLOCK, IKEEP(1,5) WILL BE SET TO -1.
C     IW    INTEGER ARRAY LENGTH 8*N.
C     .     USED AS WORKSPACE.
C     W     REAL ARRAY  LENGTH N.  USED BY F01BRQ BOTH AS WORKSPACE
C     .     AND TO RETURN GROWTH ESTIMATE IN W(1).  THE USE OF THIS
C     .     ARRAY BY M28Y12 IS THUS OPTIONAL DEPENDING ON THE VALUE
C     .     OF PARAMETER GROW.
C     LBLOCK  LOGICAL IF TRUE, F01BRY IS USED TO FIRST PERMUTE
C     .     THE MATRIX TO BLOCK LOWER TRIANGULAR FORM.
C     GROW  LOGICAL IF TRUE, THEN AN ESTIMATE OF THE INCREASE
C     .     IN SIZE OF MATRIX ELEMENTS DURING L/U DECOMPOSITION IS
C     .     GIVEN BY F01BRQ IN W(1).
C     ABORT  LOGICAL ARRAY LENGTH 4.
C     IF ABORT(1)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING STRUCTURAL SINGULARITY.
C     IF ABORT(2)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING NUMERICAL SINGULARITY.
C     IF ABORT(3)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY WHEN THE
C     .     AVAILABLE SPACE IN A/ICN IS FILLED UP BY THE PREVIOUSLY
C     .     DECOMPOSED ACTIVE, AND UNDECOMPOSED PARTS OF THE MATRIX.
C     IF ABORT(4)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING DUPLICATE ELEMENTS IN THE INPUT MATRIX.
C     IDISP  INTEGER ARRAY LENGTH 11. USED TO COMMUNICATE
C     .     INFORMATION ABOUT THE DECOMPOSITION TO THE USER AND
C     .     ALSO BETWEEN THIS CALL TO M28Y12 AND SUBSEQUENT CALLS
C     .     TO M28BSF AND F04AXF.
C     ON EXIT -
C     IDISP(1) AND IDISP(2) INDICATE THE POSITION IN ARRAYS A AND
C     .     ICN OF THE FIRST AND LAST ELEMENTS IN THE L/U
C     .     DECOMPOSITION OF THE DIAGONAL BLOCKS RESPECTIVELY.
C     IDISP(3)= IRNCP  THE NUMBER OF COMPRESSES ON ARRAY IRN.
C     IDISP(4)= ICNCP  THE NUMBER OF COMPRESSES ON ARRAYS ICN/A.
C     IDISP(5)= IRANK  ESTIMATED RANK OF THE MATRIX.
C     IDISP(6)= MINIRN MINIMUM LENGTH OF ARRAY IRN FOR SUCCESS ON
C     .     FUTURE RUNS.
C     IDISP(7)= MINICN MINIMUM LENGTH OF ARRAYS ICN/A FOR SUCCESS
C     .     ON FUTURE RUNS.
C     IDISP(8)= NUMNZ  STRUCTURAL RANK OF MATRIX.
C     IDISP(9)= NUM    NUMBER OF DIAGONAL BLOCKS.
C     IDISP(10)=LARGE  SIZE OF LARGEST DIAGONAL BLOCK.
C     IDISP(11)=Number of singular row.
C     IFAIL  INTEGER VARIABLE  USED AS ERROR FLAG BY ROUTINE.
C     THR    REAL VARIABLE  SET BY USER AS ZLATEV THRESHOLD
C
C     INTERNAL VARIABLES AND WORKSPACE USED IN  M28Y12 ARE DEFINED
C     WITHIN THE SUBROUTINE IMMEDIATELY PRIOR TO THEIR FIRST USE.
C     .. SCALAR ARGUMENTS ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='M28Y12')
      DOUBLE PRECISION U, THR
      INTEGER IFAIL, LICN, LIRN, N, NZ
      LOGICAL GROW, LBLOCK
C     .. ARRAY ARGUMENTS ..
      DOUBLE PRECISION A(LICN), W(N)
      INTEGER ICN(LICN), IDISP(11), IKEEP(N,5), IRN(LIRN), IW(N,8)
      LOGICAL ABORT(4)
C     ..
C     .. LOCAL SCALARS ..
C$P 1
      DOUBLE PRECISION THEMAX, UPRIV, ZERO
      INTEGER I1, I, IEND, II, IPRIV4, ISAVE, J1, J2, J, JAY, JJ,KNUM,
     * LENGTH, MOVE, NADV, NERR, NEWJ1, NEWPOS
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
      CHARACTER*65      REC(3)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          F01BRQ, F01BRR, M30Y12, F01BRY, F01BRZ, X04AAF,
     *                  X04ABF, X04BAF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MOD
C     .. Data statements ..
      DATA ZERO /0.0D0/
      ISAVE = IFAIL
      IFAIL = 0
C     NERR IS THE UNIT NUMBER FOR ERROR MESSAGES
C     NADV IS THE UNIT NUMBER FOR DUPLICATE ELEMENT WARNING MESSAGES
      CALL X04AAF(0, NERR)
      CALL X04ABF(0, NADV)
C     UPRIV PRIVATE COPY OF U IS USED IN CASE IT IS OUTSIDE
C     RANGE  ZERO TO ONE  AND  IS THUS ALTERED BY M30Y12.
      UPRIV = U
C     SIMPLE DATA CHECK ON INPUT VARIABLES AND ARRAY DIMENSIONS.
      IF (N.GT.0) GO TO 20
      IFAIL = 8
      GO TO 320
   20 IF (NZ.GT.0) GO TO 40
      IFAIL = 9
      GO TO 320
   40 IF (LICN.GE.NZ) GO TO 60
      IFAIL = 10
      GO TO 320
   60 IF (LIRN.GE.NZ) GO TO 80
      IFAIL = 11
      GO TO 320
C
C     DATA CHECK TO SEE IF ALL INDICES LIE BETWEEN 1 AND N.
   80 DO 100 I=1,NZ
         IF (IRN(I).GT.0 .AND. IRN(I).LE.N .AND. ICN(I).GT.0 .AND.ICN(I)
     *   .LE.N) GO TO 100
C     ** CODE FOR OUTPUT OF ERROR MESSAGE **************************
         IF (MOD(ISAVE/10,10).NE.0) THEN
            WRITE (REC,FMT=99999) I, A(I), IRN(I), ICN(I)
            CALL X04BAF(NERR,REC(1))
            CALL X04BAF(NERR,REC(2))
            CALL X04BAF(NERR,REC(3))
         END IF
C     ** END OF CODE FOR OUTPUT OF ERROR MESSAGE *******************
         IFAIL = 12
  100 CONTINUE
      IF (IFAIL.GT.0) GO TO 340
C
C     SORT ELEMENTS INTO ROW ORDER.
C
      CALL F01BRZ(N, NZ, A, ICN, IW, IRN)
      IPRIV4 = IW(1,1)
C     PART OF IKEEP IS USED HERE AS A WORK-ARRAY.  IKEEP(I,2) IS
C     THE LAST ROW TO HAVE A NON-ZERO IN COLUMN I.  IKEEP(I,3)
C     IS THE OFF-SET OF COLUMN I FROM THE START OF THE ROW.
      DO 120 I=1,N
         IKEEP(I,2) = 0
         IKEEP(I,1) = 0
  120 CONTINUE
C
C     CHECK FOR DUPLICATE ELEMENTS .. SUMMING ANY SUCH ENTRIES AND
C     PRINTING A WARNING MESSAGE ON UNIT NADV.
C     MOVE IS EQUAL TO THE NUMBER OF DUPLICATE ELEMENTS FOUND.
      MOVE = 0
C     THE LOOP ALSO CALCULATES THE LARGEST ELEMENT IN THE MATRIX,
C     THEMAX.
      THEMAX = ZERO
C     J1 IS POSITION IN ARRAYS OF FIRST NON-ZERO IN ROW.
      J1 = IPRIV4
      DO 220 I=1,N
         IF (I.NE.N) GO TO 140
         IPRIV4 = NZ + 1
         GO TO 160
  140    IPRIV4 = IW(I+1,1)
  160    LENGTH = IPRIV4 - J1
         IF (LENGTH.EQ.0) GO TO 220
         J2 = IPRIV4 - 1
         NEWJ1 = J1 - MOVE
         DO 200 JJ=J1,J2
            J = ICN(JJ)
            THEMAX = DMAX1(THEMAX,DABS(A(JJ)))
            IF (IKEEP(J,2).EQ.I) GO TO 180
C     FIRST TIME COLUMN HAS OCURRED IN CURRENT ROW.
            IKEEP(J,2) = I
            IKEEP(J,3) = JJ - MOVE - NEWJ1
            IF (MOVE.EQ.0) GO TO 200
C     SHIFT NECESSARY BECAUSE OF  PREVIOUS DUPLICATE ELEMENT.
            NEWPOS = JJ - MOVE
            A(NEWPOS) = A(JJ)
            ICN(NEWPOS) = ICN(JJ)
            GO TO 200
C     DUPLICATE ELEMENT.
  180       MOVE = MOVE + 1
            LENGTH = LENGTH - 1
            JAY = IKEEP(J,3) + NEWJ1
C     ** CODE FOR OUTPUT OF WARNING MESSAGE ************************
            IF (MOD(ISAVE/100,100).NE.0) THEN
               WRITE (REC,FMT=99998) I, J, A(JJ)
               CALL X04BAF(NADV,REC(1))
               CALL X04BAF(NADV,REC(2))
               CALL X04BAF(NADV,REC(3))
            END IF
C     ** END OF CODE FOR OUTPUT OF WARNING MESSAGE *****************
            IF (ABORT(4)) IFAIL = 13
C            write(0,*)' A(jay)=',A(jay),'  A(jj)=',A(jj)
            A(JAY) = A(JAY) + A(JJ)
            THEMAX = DMAX1(THEMAX,DABS(A(JAY)))
  200    CONTINUE
         IKEEP(I,1) = LENGTH
         J1 = IPRIV4
  220 CONTINUE
C     IF ABORT(4) IS .TRUE., DUPLICATE ELEMENTS ARE REGARDED AS
C     AN ERROR
      IF (IFAIL.GT.0) GO TO 320
C
C     KNUM IS ACTUAL NUMBER OF NON-ZEROS IN MATRIX WITH ANY
C     MULTIPLE ENTRIES COUNTED ONLY ONCE.
      KNUM = NZ - MOVE
      IF (.NOT.LBLOCK) GO TO 240
C
C     PERFORM BLOCK TRIANGULARISATION.
      IFAIL = ISAVE
      CALL F01BRY(N, ICN, A, LICN, IKEEP, IDISP, IKEEP(1,2),IKEEP(1,3),
     * IKEEP(1,5), IW(1,3), IW, ABORT(1), IFAIL)
C     ON EXIT FROM F01BRY IDISP(8)=NUMNZ, IDISP(9)=NUM,
C     IDISP(10) = LARGE.
      IF (IFAIL.NE.0) GO TO 340
      GO TO 300
C
C     BLOCK TRIANGULARIZATION NOT REQUESTED.
C     MOVE STRUCTURE TO END OF DATA ARRAYS IN PREPARATION FOR
C     M30Y12.
C     ALSO SET LENOFF(1) TO -1 AND SET PERMUTATION ARRAYS.
  240 DO 260 I=1,KNUM
         II = KNUM - I + 1
         NEWPOS = LICN - I + 1
         ICN(NEWPOS) = ICN(II)
         A(NEWPOS) = A(II)
  260 CONTINUE
      IDISP(1) = 1
      IDISP(2) = LICN - KNUM + 1
      DO 280 I=1,N
         IKEEP(I,2) = I
         IKEEP(I,3) = I
  280 CONTINUE
      IKEEP(1,5) = -1
C
C     PERFORM L/U DECOMOSITION ON DIAGONAL BLOCKS.
  300 IFAIL = ISAVE
      CALL M30Y12(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IDISP,IKEEP(1,2),
     * IKEEP(1,3), IRN, LIRN, IW(1,2), IW(1,3),IW(1,4), IW(1,5), IW(1,6)
     *, IW(1,7), IW(1,8), IW, UPRIV,ABORT, IFAIL, THR)
C
C     ON EXIT FROM M30Y12 IDISP(3) TO IDISP(7) CONTAIN
C     IRNCP, ICNCP, IRANK, MINIRN AND MINICN RESPECTIVELY.
      IDISP(6) = MAX0(IDISP(6),NZ)
      IDISP(7) = MAX0(IDISP(7),NZ)
      IF (IFAIL.GT.0) GO TO 340
C
C     REORDER OFF-DIAGONAL BLOCKS ACCORDING TO PIVOT PERMUTATION.
      I1 = IDISP(1) - 1
      IF (I1.NE.0) CALL F01BRR(N, ICN, A, I1, IKEEP(1,5),IKEEP(1,2),
     * IKEEP(1,3), IW, IRN)
C
C     OPTIONALLY CALCULATE ELEMENT GROWTH ESTIMATE.
      I1 = IDISP(1)
      IEND = LICN - I1 + 1
      IF (GROW) CALL F01BRQ(N, ICN, A(I1), IEND, IKEEP, IKEEP(1,4),W)
C     INCREMENT GROWTH ESTIMATE BY ORIGINAL MAXIMUM ELEMENT.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (IFAIL.EQ.0) GO TO 360
C
  320 CONTINUE
C     ** CODE FOR OUTPUT OF ERROR MESSAGES *************************
      IF (MOD(ISAVE/10,10).EQ.0) GO TO 340
      IF (IFAIL.EQ.-2) WRITE (REC,FMT=99992)
      IF (IFAIL.EQ.-1) WRITE (REC,FMT=99991)
      IF (IFAIL.EQ.8) WRITE (REC,FMT=99997) N
      IF (IFAIL.EQ.9) WRITE (REC,FMT=99996) NZ
      IF (IFAIL.EQ.10) WRITE (REC,FMT=99995) LICN, NZ
      IF (IFAIL.EQ.11) WRITE (REC,FMT=99994) LIRN, NZ
      IF (IFAIL.EQ.13) WRITE (REC,FMT=99993)
      IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.-2 .OR. IFAIL.EQ.13 .OR.
     *    (IFAIL.GE.8 .AND. IFAIL.LE.11)) THEN
         CALL X04BAF(NERR,REC(1))
         CALL X04BAF(NERR,REC(2))
      END IF
C     ** END OF CODE FOR OUTPUT OF ERROR MESSAGES ******************
C
  340 IFAIL = P01ABF(ISAVE,IFAIL,SRNAME,0,P01REC)
  360 RETURN
99999 FORMAT (/' ON ENTRY IRN(I) OR ICN(I) IS OUT OF RANGE - I = ',I8,
     *  /'  A(I) = ',1P,D14.6,'  IRN(I) = ',I8,'  ICN(I) = ',I8)
99998 FORMAT (/' M28Y12 FOUND DUPLICATE ELEMENT WITH INDICES ',I6,', ',
     *  I6,/'  VALUE = ',1P,D14.6)
99997 FORMAT (/' ON ENTRY N .LE. 0 , N = ',I10)
99996 FORMAT (/' ON ENTRY NZ .LE. 0 , NZ = ',I10)
99995 FORMAT (/' ON ENTRY LICN .LT. NZ , LICN = ',I8,'  NZ = ',I8)
99994 FORMAT (/' ON ENTRY LIRN .LT. NZ , LIRN = ',I8,'  NZ = ',I8)
99993 FORMAT (/' DUPLICATE ELEMENTS FOUND ON INPUT - SEE ADVISORY MESS',
     *  'AGES')
99992 FORMAT (/' MATRIX IS NUMERICALLY SINGULAR - DECOMPOSITION COMPLE',
     *  'TED')
99991 FORMAT (/' MATRIX IS STRUCTURALLY SINGULAR - DECOMPOSITION COMPL',
     *  'ETED')
      END
